package ncorps.fourier;

import ncorps.commun.ConditionInitiale;
import Jama.Matrix;

public class Interpolation {
//http://jgarrigues.perso.ec-marseille.fr/efhtml/coursefnode10.html

    private double[][][][] PolynomeInterpDuneMaille;
    private final int a = 1;
    private double[][] PotentielSurGrille;
    private double[][] gcoord;
    private double[][] GradientPotentielPourN;
//    private double[] PotentielScalairePourN;

    public Interpolation(ConditionInitiale CI, PotentielGrille potentielGrille) {
        PotentielSurGrille = potentielGrille.getPotentiel();
        gcoord = CI.getGcoord();
        GradientPotentielPourN = new double[CI.get_N()][2];
//        PotentielScalairePourN = new double[CI.get_N()];
        PolynomeInterpDuneMaille = new double[4][4][ConditionInitiale._XDimG][ConditionInitiale._YDimG];
        for (int x0 = 0; x0 < ConditionInitiale._XDimG; x0++) {
            for (int y0 = 0; y0 < ConditionInitiale._YDimG; y0++) {
                initPolynomeDuneMaille(x0, y0);
            }
        }
        double x, y;
        int x0, y0;
        for (int n = 0; n < CI.get_N(); n++) {
            x = gcoord[n][1];
            y = gcoord[n][2];
            x0 = (int) x;
            y0 = (int) y;
            gradientPotentiel(x0, y0, x, y, n);
//            Potentiel(x0, y0, x, y, n);
//            gradientPotentiel2(x0, y0, n);
        }
    }

//    private void Potentiel(int x0, int y0, double x, double y, int n) {
//        if (inBounds(x0, y0, ConditionInitiale._XDimG, ConditionInitiale._YDimG)) {
//            double potentiel0 = PotentielSurGrille[x0][y0] * EvalPolynInterpolation(x, y, x0, y0, 0, 0);
//            double potentiel1 = PotentielSurGrille[x0 + a][y0] * EvalPolynInterpolation(x, y, x0, y0, 1, 0);
//            double potentiel2 = PotentielSurGrille[x0 + a][y0 + a] * EvalPolynInterpolation(x, y, x0, y0, 2, 0);
//            double potentiel3 = PotentielSurGrille[x0][y0 + a] * EvalPolynInterpolation(x, y, x0, y0, 3, 0);
//            PotentielScalairePourN[n] = potentiel0 + potentiel1 + potentiel2 + potentiel3;
//        }
//    }

    private void gradientPotentiel2(int x0, int y0, int n) {
        if (inBounds(x0, y0, ConditionInitiale._XDimG, ConditionInitiale._YDimG)) {
            GradientPotentielPourN[n][1] = ((PotentielSurGrille[x0][y0 + 1] + PotentielSurGrille[x0 + 1][y0 + 1]) / 2) - ((PotentielSurGrille[x0][y0] + PotentielSurGrille[x0 + 1][y0]) / 2);
            GradientPotentielPourN[n][0] = ((PotentielSurGrille[x0][y0 + 1] + PotentielSurGrille[x0][y0]) / 2) - ((PotentielSurGrille[x0 + 1][y0 + 1] + PotentielSurGrille[x0 + 1][y0]) / 2);
        }
    }

    private void gradientPotentiel(int x0, int y0, double x, double y, int n) {
        if (inBounds(x0, y0, ConditionInitiale._XDimG, ConditionInitiale._YDimG)) {
            //variation selon X en X,Y du potentiel
            double potentiel_xy0 = PotentielSurGrille[x0][y0] * EvalPolynInterpolation(x0, y, x0, y0, 0, 0);
            double potentiel_xy1 = PotentielSurGrille[x0 + a][y0] * EvalPolynInterpolation(x0, y, x0, y0, 1, 0);
            double potentiel_xy2 = PotentielSurGrille[x0 + a][y0 + a] * EvalPolynInterpolation(x0, y, x0, y0, 2, 0);
            double potentiel_xy3 = PotentielSurGrille[x0][y0 + a] * EvalPolynInterpolation(x0, y, x0, y0, 3, 0);
            double potentielYA = potentiel_xy0 + potentiel_xy1 + potentiel_xy2 + potentiel_xy3;
            potentiel_xy0 = PotentielSurGrille[x0][y0] * EvalPolynInterpolation(x0 + a, y, x0, y0, 0, 0);
            potentiel_xy1 = PotentielSurGrille[x0 + a][y0 + a] * EvalPolynInterpolation(x0 + a, y, x0, y0, 1, 0);
            potentiel_xy2 = PotentielSurGrille[x0][y0 + a] * EvalPolynInterpolation(x0 + a, y, x0, y0, 2, 0);
            potentiel_xy3 = PotentielSurGrille[x0][y0 + a] * EvalPolynInterpolation(x0 + a, y, x0, y0, 3, 0);
            double potentielYB = potentiel_xy0 + potentiel_xy1 + potentiel_xy2 + potentiel_xy3;
            GradientPotentielPourN[n][0] = potentielYB - potentielYA;
            //
            //variation selon Y en X,Y du potentiel
            potentiel_xy0 = PotentielSurGrille[x0][y0] * EvalPolynInterpolation(x, y0, x0, y0, 0, 0);
            potentiel_xy1 = PotentielSurGrille[x0 + a][y0] * EvalPolynInterpolation(x, y0, x0, y0, 1, 0);
            potentiel_xy2 = PotentielSurGrille[x0 + a][y0 + a] * EvalPolynInterpolation(x, y0, x0, y0, 2, 0);
            potentiel_xy3 = PotentielSurGrille[x0][y0 + a] * EvalPolynInterpolation(x, y0, x0, y0, 3, 0);
            double potentielXA = potentiel_xy0 + potentiel_xy1 + potentiel_xy2 + potentiel_xy3;
            potentiel_xy0 = PotentielSurGrille[x0][y0] * EvalPolynInterpolation(x, y0 + a, x0, y0, 0, 0);
            potentiel_xy1 = PotentielSurGrille[x0 + a][y0 + a] * EvalPolynInterpolation(x, y0 + a, x0, y0, 1, 0);
            potentiel_xy2 = PotentielSurGrille[x0][y0 + a] * EvalPolynInterpolation(x, y0 + a, x0, y0, 2, 0);
            potentiel_xy3 = PotentielSurGrille[x0][y0 + a] * EvalPolynInterpolation(x, y0 + a, x0, y0, 3, 0);
            double potentielXB = potentiel_xy0 + potentiel_xy1 + potentiel_xy2 + potentiel_xy3;
            GradientPotentielPourN[n][1] = potentielXB - potentielXA;
        }
    }

    private void initPolynomeDuneMaille(int x0, int y0) {
        PolynomeInterpDuneMaille[0][0][x0][y0] = 1;
        PolynomeInterpDuneMaille[0][1][x0][y0] = 1;
        PolynomeInterpDuneMaille[0][2][x0][y0] = 1;
        PolynomeInterpDuneMaille[0][3][x0][y0] = 1;
        //
        PolynomeInterpDuneMaille[1][0][x0][y0] = x0;
        PolynomeInterpDuneMaille[1][1][x0][y0] = x0 + a;
        PolynomeInterpDuneMaille[1][2][x0][y0] = x0 + a;
        PolynomeInterpDuneMaille[1][3][x0][y0] = x0;
        //
        PolynomeInterpDuneMaille[2][0][x0][y0] = y0;
        PolynomeInterpDuneMaille[2][1][x0][y0] = y0;
        PolynomeInterpDuneMaille[2][2][x0][y0] = y0 + a;
        PolynomeInterpDuneMaille[2][3][x0][y0] = y0 + a;
        //
        PolynomeInterpDuneMaille[3][0][x0][y0] = x0 * y0;
        PolynomeInterpDuneMaille[3][1][x0][y0] = (x0 + a) * y0;
        PolynomeInterpDuneMaille[3][2][x0][y0] = (x0 + a) * (y0 + a);
        PolynomeInterpDuneMaille[3][3][x0][y0] = x0 * (y0 + a);
    }

    private double EvalPolynInterpolation(double X, double Y, int x0, int y0, int I, int J) {
        double[][] PMailleEnCours = new double[4][4];
        for (int i = 0; i < 4; i++) {
            for (int j = 0; j < 4; j++) {
                PMailleEnCours[i][j] = PolynomeInterpDuneMaille[i][j][x0][y0];
            }
        }
        Matrix poly = new Matrix(PMailleEnCours);
        poly = poly.inverse();
        double[] TVectBase = new double[4];
        TVectBase[0] = 1;
        TVectBase[1] = X;
        TVectBase[2] = Y;
        TVectBase[3] = X * Y;
        Matrix MVectBase = new Matrix(TVectBase, 4);
        Matrix valinterp = poly.times(MVectBase);
        return valinterp.get(I, J);
    }

    private boolean inBounds(int x0, int y0, int _XDimG, int _YDimG) {
        return ((x0 > 0) && (y0 > 0) && (x0 < _XDimG - 1) && (y0 < _YDimG - 1));
    }

    public double[][] getGradientPotentielPourN() {
        return GradientPotentielPourN;
    }

//    public double[] getPotentielScalairePourN() {
//        return PotentielScalairePourN;
//    }
}








